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AND  ITS  COMPUTATION  I 
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Benren  Zhu 

Coapster  Science  Department 
Yale  University 

t  Abstract 

C.. 

A 

In  tbit  paper  we  dhall  investigate,  a  type  of  nonlinear  parabolic 
equation,  coupled  equations  and  their  positive  solutions.  The  properties  of 
maximum  and  the  nonlinear  instabilities  will  be  discussed  briefly;  for 
obtaining  its  numerical  solution  a  special  scheme  wiU  be  presented*. x 

1*  Introduction 

In  the  present  paper  we  shall  investigate  a  type  of  nonlinear  parabolic 
equation  and  coupled  equations  to  which  the  similar  one  may  come  from  physics, 
chemistry  and  ecology*  Specifically,  the  coupled  equations  are  from  solid 
electronics  and  have  been  simplified  to  describe  the  dynamic  distributions  of 
two  kinds  of  charge  carriers  so  called  carrier  equations*  The  unknown 
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functions  in  these  problem  represent  the  densities  of  certsin  particles  or 
bodies.  The  positivity  of  those  functions  in  their  mathematical  nodel  will 
play  essential  role.  Ignorance  of  this  property  often  cause  serious 
instability  especially  in  nonlinear  problems,  which  we  call  nonlinear 
instability,  A  mathematical  model  describing  a  natural  phenomena  which  is 
reproducible  must  be  properly  defined  and  continuously  dependent  on  a  number 
of  given  conditions.  We  are  concerned  with  the  case  in  which  there  is  a  steady 
state  solution  and  a  transitional  solution  wherein  the  former  is  the  limit  of 
the  latter  as  time  increase  to  infinity. 

In  order  to  solve  such  a  problem  we  suggest  using  a  special  scheme 
, positive  scheme  ,  which  will  be  used  in  computation  of  the  positive  solution. 
This  assures  us  to  avoid  the  nonlinear  instability  • 

Many  authors  have  studied  this  type  of  problem.  Some  of  them  refered  to 
Lyapunov  stability,  some  of  them  refered  to  the  direct  integration  methods. 
But  we  prefer  the  positivity  of  the  solution,  and  emphasis  is  placed  on  the 
construction  of  the  computational  scheme, 

1.  A  Nonlinear  Parabolic  Equation 

Let  us  begin  with  a  single  nonlinear  parabolic  equation  as  follows: 

§2-.0-b<u-i<*))ir  (A) 

,.b  >  0.  -1  <  *  <  1.  t  >  0. 

whet,  D(x,O)a0(x)i  U(+l,  t)  are  given  positive  functions  and  K(x)  is  given  as 
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an  arbitrary  continuous  function*  A  set  of  functions  which  consists  of  all 
positive  continuous  functions  we  call  C+ ;  and  we  call  such  a  solution  a 
positive  solution. 


THEOREM  1*  If  U(x,t)  is  a  positive  solution  of  (A)  then  the  following 
inequality  holds: 

0  <  U(x,t>  i  I  -  Max  (MaxK(x) 0  Max  U}  (I) 

X  OT^tl 

Proof.  On  the  contrary  there  Boat  be  x  e  (-1.1)  and  t  >0  anch  that 


Max.  U(x.t)  >  X.  ft>  1  0.  0. 

tit  9t  dx 

and  (U(x  ,t)-K(x*))B(x*,t)  >  0.  It  confiiota  with  the  equation  (A). 


THEOREM  2.  The  positive  aolntion  of  (A)  ia  unique. 


The  proof  of  theorem  2  ia  oaitted.  Considering  the  ateady  atate  aolntion 
of  (A)  it  ia  obvions  that  thia  aolntion  ia  free  of  variable  t  and  ia  the  limit 
of  the  aolntion  V(x,t)  daring  t — >*».  The  ateady  atate  aolntion  of  (A) 
aatiafiea  the  nonlinear  elliptic  equation  and  boundary  conditions  which  we 
call  problem  (A*).  First  we  divide  the  function  K(x)  into  two  positive 
functions  K*(x)  and  I  (x)  anch  that  K(xX*(x)-lT(x)  i.e. 


r+.  .  .  K<x)  +  p(x)  ...  if  X  >  0; 

1  p(x)  .  otherwise; 


X'tx) 


p(x)  .  if  X  >  0; 

**K(x)tp(x)  ...  otherwise. 


where  p(x)  ia  an  arbitrary  non-negative  function. 
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Next  we  suggest  tn  iterative  process  which  is  very  important  in  finding  a 
solution  of  (A).  Consider  the  following  iterative  process: 

a|~-b  [  C,1DB"1-K+U1'"1  +K  V] 

dt  ®x  n-1,2 .  (A.) 

where  D*(+l,t)«U(+l,t)  ,  D"(x,0)«l>(x) ,  but  U°  is  any  function  satisfying  (I)> 
Each  step  of  the  prooess  is  forwarded  by  solving  a  linear  parabolic  equation. 

THEOREM  3.  All  iterative  solutions  of  (An)  satisfy  the  following 

inequality 

0  <  D*(x,t)  IK*  Max[MaxK+(x) ,  Max  0] 

°  x  or^Stl 

Proof.  For  n-0  theorea  is  true  .  We  suppose  inequality  valid  for  all 

integers  less  than  n  .  Then  we  check  the  inequaliry  for  integer  n.  On  the 

$  $  ■  t  i 

contrary  there  east  be  (z  »t  )  in  the  region  such  that  function  u  (x  #t)  at  a 
aaxiaua.  Both  of  then  conflict  with  equality  of  (A^) . 

Now  we  discuss  the  convergence  of  the  sequence  (1)*)  in  the  functional 
space  L*.  Set  s"  ,  then  a*' a  satiafy  the  following  equation  with 

hoaogeneoua  initial  and  boundary  conditions: 

Multiply  both  sides  by  eB  and  integrate  on  whole  interval  (-1,1).  using 
the  following  notations: 

|v(x)l*[f  IvMl^dx]1^2; 

•'-1 

r*^ 

(u,v)-  I  u  V  dx. 

-1 


A 
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Then  we  have 


<m>|i|.-i2-.ii£.*i2-»i,  «r2.*i2- 

-bl/  K"eB|2+b((K+-DB"1)e“_1,eB)i 


i  -a|eB|+bC|eB|  le*-1 1 . 


By  asing  the  well  known  inequality  lei  £  Ir^el  we  have 

ox 


where  n1“2a-bC,  n*®bC,  i.e. 

3t 


i-»tU“12+o»t."_1l2 


s|(l.*l2.***)  i  .*l.*-1l2.n*t. 


thus  we  have 


i.*i2.*1*  i  «*f  i.*-ii2.**tit  x 

J0 

i  (».)2jVi.-2i2."‘''d.'d.  i ... 

i  l]f,(t-.l*"2l.1l2."**dt  i 

*  A 


i  [(n«t)B"1/(«-l)!]  Mar  le^V*1) 

Therefore  the  following  theorem  oan  be  drived. 

Theorem  4.  Iterative  aeqnenoe  (D1)  is  convergent  in  any  finite  region  in 
the  sense  of  the  L*-norm. 


Further  results  about  the  generalised  eolation  of  (A)  are  beyond  oar  main 
purpose.  Vs  are  interested  in  this  iterative  proeess  (A  )  and  its  motivated 
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scheae  of  computation. 


2.  Nonlinear  Instability  and  A  Computational  Scheme 

When  the  solutions  of  (A)  are  confined  in  C+,  then  the  positive  solution 
of  (A*)  could  be  obtained  as  a  limit  of  (A)  during  t — >®.  We  call  stable 
steady  state  solution.  But  once  ve  disregard  this  restriction,  the  unstable 
phenomena  will  often  occur.  For  example  ,  if  EsO*const  when  <  0;  a=0,  then 
as  the  explicit  solution  shows: 

U»C/[1-{[(U  -C)/U  ]e"bCt}]— >•  as  t—  >[l/(bC)]Ln[(U  -C)/UJ. 

o  o  o  o 


The  instability  becomes  obvious.  According  to  variety  of  the  constants  C 
and  the  possibilities  are  listed  in  the  table  1  below  and  show  the  stable 
(unstable)  region  in  figure  1. 


I  C  <  0  I  c  >  o  I 
I - 1 - 1- 


I  stable  I  stable  I  0  >  0  I 

| - 1 - 1 - 1 

I  Halted  | unstable  |C<  U  <  Ol 

| - 1 - 1 - 1 

lonstablelonstablel  D  <  C,0| 


stable 


stable 
Figure  1. 


Table  1. 

The  very  unfavorably  unstable  situation  will  develop  if  one  puts  the  initial 
value  in  the  unstable  region.  So  ve  aust  be  careful  to  construct  the 
approxiaate  soheae  so  that  once  any  initial  function  is  selected  in  the  stable 
region  then  all  the  interaidate  results  will  never  be  out  of  the  stable 
region.  Thus  the  coaputation  can  be  possibly  oarried  on  until  a  solution  has 
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been  found* 


Now  wo  cone  to  a  computational  scheme  for  problem  (A),  First  of  all  we 
divide  in  terval  (-1,1)  into  2N  subintervals  by  nodes: 

-1^_N<*-N+1<  •  •  •*-2<*-1  <vv  •  •  *J Vi<: V1  • 

Let  h^-x^-x^gi—N+l#.  ..,-1,0,1, .  ..,N  be  lengths  of  these  subintervals  and 
use  the  notations: 


u“«U(xJ#t  ),  t  smAt,  bf0,1,2,... 
i  l  m  m 


1  2h.  2h. 

- 1 — 

ki Vi  Wi  Wi 


where  At  is  the  length  of  the  subinterval  in  t-axis. 


We  set  the  difference  equation  as  follows: 


(Ah) 


where  oJ|,D^,U°  *re  positive,  i**0,±l,. .  ,+(N-l) ,  w»l,2,... 

Let  Ch+  be  e  eet  of  positive  functions  defined  on  these  nodes.  We  call  a 
solution  of  (Ah) ,  which  belong  to  Ch+,  the  positive  solution. 

Theorea  5.  The  coefficient  aatriz  of  the  unknown  functions  in  (Ah)  has 
a  unique  inverse  with  non-negative  eleaents. 

It  is  easy  to  see,  if  we  reaove  all  the  teras  of  to  the  left  sides  of 
the  equations  and  O^^/At  to  the  right  sides,  then  the  coefficients  of  D^J  fora 
a  speoial  aatriz  with  doainate  diagonal  eleaents  negative,  so  there  is  a 


nonnegative  inverse* 


Theorem  6*  Under  the  conditions  of  theorem  3  the  equation  (Ah)  has  the 
solution  in  Ch+  and  the  following  inequality  holds: 

o  <  ^  X  <•  *M  ”!>• 

i  or  I*+N 

The  proof  is  much  the  same  as  for  section  1* 


Theorem  7*  If  the  solution  of  (A)  has  all  continuous  derivatives  for  x 
and  t  required,  then  the  difference  solution  approximates  differential 
solutions  in  any  finite  section  of  the  region  in  the  sense  of  norm* 

In  fact*  suppose  we  take  e*  as  the  difference  between  the  solution  of  (A) 
and  the  solution  of  (Ah).  Then  e*  suits  the  following  difference  equations: 

(e®-c^~*)/At«a6^e“+ 

+bK*e^“1-bK~e*-bU^~1 e“-bu“e^~1+R^ 

R®- approximate  error*  Multiply  the  two  aides  of  the  equation  by  e*  h^,  where 
h^*(h^+h^+^)/2  and  sum  up  from  i=-N+l  to  N-l.  We  introduce  the  notations  about 
mesh  functions  in  the  same  way  as  the  continuous  functions 


then  we  obtain 


l(^I|^1V(vr*?»/V2- 

le*lJ+«At  U“|J+bAt I (l“)1/2e"|»- 

!l  1 •T"1  { i+bAt  (i+_d*)  1 + 
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Other  schemes  may  be  available »  but  scheme  (Ah)  is  special  from  the 
iterative  process  (Am)  in  section  1.  Under  a  slightly  strict  condition  we 
obtain  the  stability  of  the  difference  equations  (Ah)  with  regard  to  initial 
values. 

Theorem  8.  If  a  >  b(Ma x  K+ ) ,  then  the  difference  equations  are  stable 
with  regard  to  initial  values. 

This  means  if  U^  and  V^are  two  solutions  of  (Ah)  with  diff  different 
initial  values  ,  then  the  difference  between  them  does  not  increase.  The 
proof  is  the  same  as  for  theorem  7  when  omitting  R®. 

3.  A  Coupled  Nonlinear  Parabolic  System 


C+  we  still  ase  the  functions  K+  and  K  in  section  1.  and  denote  the  inner 


maximum  valoe  as  follows: 


U  .  *U(x  ,t  )*  Max{Max  U#Max  V] 
max-in 

*  *x  1 

-1  <  x  <  1,  t  >0 


the  same  to  V 


max- in* 


Theorem  9,  Only  in  region  K(x)  }_  0  U(x,t)  coaid  reach  U  .  and  for 

max— in 

which  the  following  inequality  holds: 

0i  U  .  -V(x,t)  i  Kp  =  MaxK+(x);  (II) 

max- in 

x 

only  in  region  K(x)  <  0  V(x,t)  coaid  reach  V  .  and  for  which  the  following 

max— in 

inequality  holds: 


0£  V  J  -U(x,t)  i  Km  *  MaxK~(x), 
max-in 

x 

where  the  functions  in  each  inequality  are  evaluated  at  the  same  point, 


(II) 


(Bm) 


Using  the  same  idea  as  in  section  1  we  suggest  an  iterative  process: 

if**  ■ 

ifv"  ■  *25?v*-  V  '1-  yV-f-riv--1 

where  U^V®  mF=0,l,2,...  have  the  same  initial-boundary  values  as  U,V  and  U°»V° 
are  a  pair  of  arbitrary  functions  satisfying  (II)  in  C+. 


Theorem  10#  Each  pair  of  the  functions  in  sequence  {U®,V®}  generated  by 
iterative  process  (Bm)  remain  in  C+  and  the  following  inequalities  hold: 

0  i  C^i.-v*  i  **» 


(Ilm) 


We  leave  the  proofs  of  these  two  theorems  to  the  readers.  Even  more  resalts 


about  convergence  and  generalized  solutions  could  be  obtained,  But  we  are 
going  to  describe  the  problem  of  instability  and  a  computational  scheme  in  the 
next  section. 


Figure  2. 

4.  Nonlinear  Instability  And  Computational  Scheme 

First  we  investigate  the  instability  of  system  (B) .  Assume  a^  =  0,  B^l , 
i*l,2,  K(x)°  const  divide  the  first  equation  by  the  other  ,  then  we  obtain  an 
integration  of  U*V  -const  where  the  constant  is  decided  by  the  initial  value 
of  U  and  V.  All  the  possibilities  about  stability  are  shown  in  figure 
2,(a)»(b),  (c).  The  arrows  represent  the  behavior  of  the  solutions  U,V  during 
t — >®,  Only  those  points  which  have  two  arrows  on  both  sides  pointing  to  them 
represent  the  stable  steady  state  solutions.  For  our  purposes  only  the 
positive  solutions  are  available  ,  i.e.  U»V  >  0,  otherwise  nonlinear 
instability  may  arise  as  seriously  as  in  section  2. 

Our  next  task  is  to  give  the  computational  scheme  of  equations  (B) .  Using 
the  same  division  of  time-space  coordinates  and  the  same  notations  as  in 
section  2,  we  suggest  the  differenoe  equations  as  follows: 
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*  *“1/  *  ®\  .v  ir+  *“1 


(vTr1)/At“#26Jvb2KiV 

-  m-1 ,  d  mx  ..  -  m~l 

"b2vl  (V°i)b2KiVi  • 


(Bh) 


where  u®^,  v®^;  u°*  v°*  i»0#+l,+2,..  are  given,  positive  and  satisfying  (Ilh) . 
We  call  a  solntion  of  (Bh)  which  belong  to  Ch+  positive  eolation. 


We  denote 

,  «u®*«  ***(  Max  Max 
max- in  i  i  ^  i  1 

v®  =v"***  Max(  Max  a".  Max  v®) 
max- in  i  .  i  .  i 

e  ee1  1 

-N  <  i  ,  i  <  +N, 

as  inner  maxima**  of  the  eolations. 

Theorem  11.  Problem  (Bh)  has  a  unique  positive  solution  in  Ch+. 

m  m 

Proof.  Remove  all  the  terms  containing  u^  or  v^  to  the  left  sides  of  the 

equations  and  u^~* ,  v^~*  to  the  right  sides.  Then  all  the  terms  on  left  sides 

mm 

are  non-negative.  The  coefficient  matrix  of  unknown  functions  o^,  v^  has 
dominant  diagonal  elements  and  negative  non-xero-off-diagonal  elements*  so  it 
has  an  unique  non-negative  inverse.  Thus  we  get  a®,  v®  £  0.  Further  we  shall 
verify  a®,  v*  >  0.  If  for  some  (m*i)  u*  or  v*  *  0*  suppose  u^^O  and  VV-  >0 
for  all  ■  <  a, ,  then  cheek  the  signs  of  every  tern  in  (Bh) ,  we  shall  see  the 
right  side  is  positive  which  is  conflict  with  left  side. 

Theorem  12.  Only  in  region  where  K(x)  1  0  oonld  u*  reach  «* 


and  to 
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which  wc  have 


Oil  .  -v"  i  I  ; 
sax- in  i  i’ 


(Ilh) 


only  in  region  where  K(x)  i.  0  could  v.  reach  v  .  and  we  have 

i  aax-in 


0  *  wl  *  V 


(Ilh) 


The  proof  of  theorem  12  can  be  carried  on  in  the  same  way  as  theorem  3  or 


10. 


Many  other  schemes  have  been  found.  Here  we  only  mention  the 
computational  scheme  (Bh) .  The  reason  is  that  does  this  scheme  (Bh)  keep 
almost  same  characteristics  as  the  original  problem  (B) .  Besides  the  equations 
of  problem  (B)  contain  the  main  core  of  many  more  complex  problems.  We  hope 
that  it  will  be  useful  to  those  who  are  interested  in  fields  we  mentioned,  and 
the  method  of  constructing  (Bh)  could  be  easily  used  in  solving  some  other 
related  problems. 


5.  Discussion  of  The  Algorithm 

We  write  equations  (Bh)  in  a  matrix  form  as  follows: 


<5^  -  S^/At 


•l 0 

Vi  » 

0 

x  i 

lo  !>,*! 

n- 


.  m-1  .  *-l 

Vi  -Vi 

vt° 

v  -  m-1 

Vi 

0  vl 

n- 


where  uj  i.  a  colain  veotor  (a*,v*)*.  Using  notations: 
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prV<vki*i)'  qi*1_pi' 


A  J'aiVi  0  r  J~Wi  0 
1  °  ”*2qlri  i,0"Viri 


B" 


F"= 

i 


l+2«1ri+b1K~At+b1n^"1At  -bjO^At 
-bj^A*  l+2VlrtjI*At*bj.flAt 
(l+bjKjAtJu^-1 

(Irt^At)*^1 


we  get  the  following  block-tridiagonal  system: 

*  *??  *  ci^«  - 
i*0,±l,±2,#..,£(N-l) 

where  are  given.  The  direct  method  can  be  used  in  solving 

eliminating  we  transform  (III)  into  (IV): 


*T  *  <5?«  ■ 


i-0,±l.+2,...,+(N-l), 


Si+l  "  <B-N+1)  *Sj+l’  ^-N+l*  *  ®-N+l  *  #FVi; 


Ci+1  "  <Bl+rWi>  Ci+1* 


*  *Bi+l“Ai+lCi*  ^Bi+l"Ai+lFi  * 
i— N+l , . . . ,-2 .-1 ,0.1 ,2 , . . . ,N-1 . 


(Ill) 
(III).  By 


(IV) 


V*  recall  that  proceaaea  (IV)  and  (V)  are  coapntable  and  the  nmber  of 
operations  is  about  0(N)  for  each  step. 
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fe  have  worked  on  this  coupled  equations  for  a  while.  Notice  1)  Theoren 
12  is  not  sufficient  for  the  stability;  2)  Theoren  12  only  provide  a  reasonable 
restriction  on  the  difference  solution  as  theoren  10  does;  3)because  of 
nonlinear  effect  At  should  not  be  too  big#  so  that  a  stable  approxination  can 
be  used. 

References 

[llNonl inear  Problens  in  the  Physical  Sciences  and  Biology 
Lecture  Notes  in  Mathenatics  (1972)  Springer-Verlag 

[2]R.  0.  Richtnyer  and  K.  V.  Norton 
Difference  nethods  for  initial  value  problens 
New  Tork#Wile3r-Interscience(1957) 


A 


